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ABSTRACT 

In previous papers we discussed results from fully time-dependent radiative transfer mod- 
els for core-collapse supernova (SN) ejecta, including the Type Il-peculiar SN 1987A, the 
more "generic" SN Il-Plateau, and more recently Type Ilb/Ib/Ic SNe. Here we describe the 
modifications to our radiative modeling code, CMFGEN, which allowed those studies to be un- 
dertaken. The changes allow for time-dependent radiative transfer of SN ejecta in homologous 
expansion. In the modeling we treat the entire SN ejecta, from the innermost layer that does 
not fall back on the compact remnant out to the progenitor surface layers. From our non-LTE 
time-dependent line-blanketed synthetic spectra, we compute the bolometric and multi-band 
light curves: light curves and spectra are thus calculated simultaneously using the same phys- 
ical processes and numerics. These upgrades, in conjunction with our previous modifications 
which allow the solution of the time dependent rate equations, will improve the modeling of 
SN spectra and light curves, and hence facilitate new insights into SN ejecta properties, the SN 
progenitors and the explosion mechanism(s). CMFGEN can now be applied to the modeling of 
all SN types. 
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1 INTRODUCTION 

Supernova (SN) spectra potentially contain a great deal of informa- 
tion about the progenitor star, about the explosion dynamics, and 
nucleosynthesis yields. Extracting this information, particularly at 
early times, is difficult due to their low densities and high veloc- 
ities at their effective photosphere. Because of the low densities, 
radiative processes tend to dominate over collisional processes and 
hence local-thermodynamic equilibrium (LTE) cannot be assumed. 
Instead we must solve the statistical equilibrium equations, and this 
requires vast amounts of atomic data, much of which has only be- 
come available over the last decade. The high velocities blend spec- 
tral features together, often making line identifications difficult, and 
this hinders spectral analysis since accurate spectral modeling is 
needed in order to interpret weak, but important, diagnostics. 

An inherent feature of SNe is that their spectra evolve with 
time, on a timescale comparable to, or shorter than, the age of 
the SN. Thus a crucial question is whether this time dependence 
needs to be taken into account when modeling SN spectraQ Can 
reliable results from spectral fitting be obtained by modeling only 
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^ When modehng light curves there are no controversies — time depen- 
dency is essential for modeling the light curves of SNe. An excellent dis- 
cussion of some of the issues regarding time-dependence is provided by 
iPinto & Eastmaij )2000l) . A related issue i s the necessity of includin g all 
terms of order v/c 'm the transfer equation; iMihalas & Mihala^ jl984l) pro- 



the photospheric layers, in the spirit of radiative-transfer studies 
on stellar atmospheres? That is, can we ignore time-dependent 
terms, and can we impose a lower boundary condition and adjust 
this boundary condition, and photospheric abundances, until our 
model spectrum matches what is observed? This approach has been 
adopted by many w o rkers i n the fi eld, following the first pa pers by 



[Branch et alj ( ll985h,lLucvl (Il987l) aiidlMazzah et al.l ( 1 19931) . 

For example, iDessart & Hillietl €00^ used the technique to 
explore the expanding photosphere method (EPM) for finding dis- 
tances to Type II SNe. With these investigations they were able to 
investigate the physics influencing the EPM technique, and hence 
understand why it was difficult to ac curately estimate t he co rrec- 
tion factors using photometry alone. iDessart & Hillietl fcOOfih ap- 
plied the EPM technique to SN 1999em, and showed that the dis- 
tance obtained agre ed within 10% of the cepheid distance. Later, 
IDessart et al.l ( 1200 8l) successfully applied the EPM method to de- 
termine the distance to SN 2005cs and SN 2006bp. Our modeling 
of these three Type II-P SNe did, however, reveal a problem that is 
also seen in modeling of other hydrogen-rich SNe — the inability 
to predict Ha during the recombination epoch (after approximately 
day 20 in Type II-P SNe). At this and later times the model Ha line 
is consistently narrower and weaker than observed. 

The Ha problem is also seen for SN 1987 A; however in 



vide an in-depth discussion of the impoitance of various terms in the trans- 
fer equation. 
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this case problems in fitting Ha occur as early as day 3 (e.g., 
ISchmutz et al. due to the more rapid coolin g of th e ejecta 

from a more compact progenitor ISchmutz et al ] (Il990h postu- 
lated that clumping might be the source of the Ha discrepancy. 
lEastman & KirshneJ ( Il989h have also modeled SN 1987 A (up to 
day 10) and did not find any Ha dis crepancy. Two effect s migh t 
contribute to this difference. First, lEastman & Kirshtieil d 19891) 
treated many metals in LTE wh ich could affect the H ioniz ation 
balance. Second, as shown by lEastman & Kirshne3 l ll989l) . the 
Ha emission strength is sensitive to the density exponent n de- 
fined by 71 = —d\np/d\nr. Low n (i.e., n — 7) give much 
strong er Ha emission than higher n (i.e., n = 20). ISchmutz et al] 
id) used a powe r law with a density exponent of 10, while 
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lEastman & Kirshnej ([19891) used a density exponent of 9. In the 
study of iDessart & HillieJ ilOOdt) it was found that while lowering 
the density exponent down from 10 to 8, and even to 6, led to an 
increase in the Ha strength it also produced a very severe mismatch 
of the model with observations of the Ca II line s. 

A related problem is that ISchmutz et al.l ( Il990t) were unable 
to fit the He I emission in SN 1987A - models predicted no He I 
A5876 although it is easily seen in the observations (up to day 4). 
Simple model changes did not improve the predictions, and hence 
they postulated an extra s ource of UV photons. The models of 
lEastman & KirshneJ ( fT989h also underestimated He I A 5876. They 
suggested a combination of effects as the likely cause — high he- 
lium abundance. X-rays produced in an interaction zone, density 
effe cts, and time-depende nce effects. 

iMitchell et"al] hoOll) found that models for 1987 A computed 
using the code PHOENIX jHauschildt & BaronI 1 19991 : IShort et al] 
Il999h . also predicted Ha emission much weaker than that ob- 
served. To explain these discrepancies they invoked mixing of 
^^Ni into the SN's outer layers — mixing at velocities (v > 
5000 km s^'^) l arger than predicted in curren t hydrodyna mic mod- 
eling jKifonidi s et al. 2006; Jogge rsTet alj [2009: Ham mer et al.l 
I2OI0IFI . The Hi and He I lines would be excited by non thermal 
electrons produced by the degr adation of gamma- rays created by 
the decay of radioactive nickel. iBaron et"al] ( |2003|) also used mix- 
ing of ^^Ni into the SN's outer layers to explain the H I and He I 
lines in early-time spectra of SN 1993W (Type II-P). 

The H I and He I discrepancies are important — either there 
is a fundamental problem with our current spectral modeling of 
SNe, or there is crucial SN physics missing from current hydro- 
dynamic models which leads, for example, to an underprediction 
of the amount and extent of '"''Ni mixing. In order to distinguish 
these possibilities it is important that we relax assumptions, s uch as 
steady state, in our SN models. Indeed, lUtrobin & Chugail ( |2005|) 
showed that it was important to include time dependence in the rate 
equations for SN 1987A. When such terms are included, it is possi- 
ble to explain the large strength of the B aimer lines without invok- 
ing non-thermal excitation resulting from the mixing of radioactive 
nickel from below the hydrogen-rich envelope. This conclusion has 
been confirmed bv lOessart & Hillieil ( l2008h who showed that inclu- 
sion of the time dependent terms allowed the Ha line strength to be 
better matched in SN 1999em, particularly at late times (>20 days). 

The necessity of including time-dependent terms is unfortu- 



^ The recent 3D simulations of iHammer et al.l )2010l) do produce nickel 
mixing out to 4500 km s~^, but this is insufficient. Further, in order to ex- 
plain Ha, which is the most significant spectral feature in Type II SN spec- 
tra, isolated bullets will not suffice - there has to be global mixing into the 
hydrogen photosphere. 



nate. With steady-state models we can readily adjust the luminos- 
ity, abundances and the density structure to tailor match the ob- 
servations. Instead we have an initial value problem requiring us 
to begin the calculations at a sufficiently early time such that the 
adopted initial model allows the full time evolution to be modeled. 
On the other hand, this necessity will allow us to better test both 
hydrodynamical models of the explosions, and models for the pro- 
genitor evolution. 

More recently iDessart & Hillied ( |20IO|) . extended their work 
to allow for time dependence of the radiation field. This removes 
issues with the inner boundary condition, and controversies on 
which terms can be omitted from the rad iative-transfer equation 
tearon et alJll996l : |Pinto & Eastmanll2000l) . With their fully time- 
dependent approach, they were able to successfully model the early 
evolution of SN 1987A, from <1 d up to 20 d. Overall, predicted 
spectra were in good agreement with observation, and showed the 
expected trend from a very high effective temperature at early times 
(Tgg- > 30, OOOK at 0.3 d) to a relatively constant effective temper- 
ature at later times (T^ff ~ 6000 K at 5 d). The latter is a con- 
sequence of the hydrogen ionization front which sets the location 
of the SN photosphere. One deficiency in the models was an ex- 
cess flux in the U band. This excess flux could arise from a de- 
ficiency in the model calculations (e.g., insufficient line blanket- 
ing) or from deficiencies in the initial model (e.g., wrong sized 
progenitor). This work, with time dependent radiative transfer and 
time dependent rate equations, showed that standard models for 
1987A, without ^"^Ni mixing beyond 4000 km s ^ could explain 
both the Hi and He I lines in SN 1987A. The result is in agree- 
ment with explosion models which require mixing out to ~ 3000 
to 4000 km s ~^ in order to explain the observed light curve of SN 
1987A (e.g.. IShigevama & Nomot3ll990l: iBlinnikov et al.ll2000l) . 
the appearance of hard X -ravs (e.g..lKumagai et al .1 19891) . and later 
nebular diagnostics (e.g.. lHachinger et al.l2009t) . 

Time de pendent radiatio n transfer has also been included 
in PHOENIX Ijack et al.l |2009|). as has th e solution of the time- 
dependent rate equations ( iDe et al. distinction between 
PHOENIX and CMFGEN is that PHOENIX works primarily with the 
intensity, while CMFGEN solves for the radiation field using mo- 
ments. The latter method has two advantages: It explicitly allows 
for electron scattering, and in the time-dependent approach a fac- 
tor of A^angio(~ND) Icss memory is required to store the time- 
dependent radiation field. A disadvantage is the generalization to 
multi-dimensions, and stability issues. 

In our work, we solve the radiative- transfer equation, or its 
moments, numerically. An alternative technique is to use a Monte- 
Carlo approach. A major proponent of the Monte-Carlo technique 
is Lu cy, who has d evoted considerabl e effort to its development 
(e.g.. lLucvl[T999allbl I2OO2I I2003L |2005|) . Techniques have been de- 
veloped th at allow the temperature structure of th e envelope to be 
solved for ( lLucvlll999al : lBiorkman & Woodll200ll) . Most, if not all, 
the Monte-Carlo codes use the Sobolev approximation for the line 
transfer: this sho uld be an excellen t approximation in SN atmo- 
spheres although lBaron et"aL I (fl99^ argue that the approximation 
is suspect in SNe because of the larger number of lines whose in- 
trinsic profiles overlap. While it is possible to undertake full non- 
LTE calculations with Monte-Carlo codes, such calculations are 
very time consuming. Thus it is usual to make approximations to 
determine the ioni zation structure [ e.g., L TE or a modified nebu- 
lar approximation: lAbbott & Lucvl ( 1 1 9851) 1 and/or the line source 
function (e.g., LTE). 

Early applications of Monte- Carlo methods to SNe are dis- 
cussed by iMazzali & Lucvl ( 1 19931) and were applied to the mod- 
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eling of SN la spectra by iMazzali et al. I l ll993l) . Since that time 
the code has undergone numerous extensions (e.g.. lMazzalill2000h 
and has been appli e d to i nterpreting spectra of many SNe (e.g., 
iMazzali et al1l2006l , I2OO9I) . The code is not time dependent, and 
generally uses a gray photosphere approximation for the inner 
boundary condition. Deviations in the level populations and ioniza- 
tion structure a re taken into account using a nebular approximation 
jMazzali et ai][200^) . These Monte-Carlo codes are neither LTE, 
nor fully non-LTE. Since line photons are allowed to scatter, non- 
LTE is partially treated even when LTE populations are assumed. 

A more recen t example of a Monte-Carlo code is SEDONA 
jKasen et alj200d') .SEDONA calculates emergent spectra, as a func- 
tion of time, assuming LTE level populations. Line transfer for the 
"strongest line s" (up to 0.5 m il lion) are treate d using the Sobolev 
approximation I S obolevi ( ll96(ll) : ICastoil | [l970l) while all remaining 
lines are treated using an expa nsion opacity formalism dKarp et al.l 
I1977I : lEastman & Pintdll993h . In SEDONA photons can be either 
absorbed by a line, scatter, or fluoresce (cause an emission in an- 
other transition but which has the same upper level). In the study 
of the width-luminosity relation for Type la SNe redistribution pro- 
cesses were taken into account using a two-le vel atom approach 
with a constant redistribution probability of 0.8 jKasen & WooslevI 
2007h whereas in ano ther study it was set to 1 (i.e., pure absorption; 



Wooslev et~al] |200"^ FI One advantage of the Monte-Carlo tech- 



nique is_the relative ease with which time-dependence is treated 
( ILucvI I2OO5I) . More rec ent work using Mo nte- Carlo techniques 
has be en undertaken by iMaurer et al.l ( 201 ih and Ijerkstrand et al.l 

l|) who solve the non-LTE time-dependent rate equations, and 
solve for the radiation field during the nebular phase. The approach 
allows for a more accurate treatment of the emergent flux for quan- 
titative analyses especially since, in the nebular phase, the coupling 
between the gas and the radiation is relatively weak. Monte-Carlo 
codes tend to be inefficient in ID but have the advantage that they 
can be readily extended to 2D and 3D. Given the complexities of 
SN modeling it is essential that we have alternative approaches 
whose results can be compared. 

In the present paper we outline the assumptions underlying our 
time-dependent radiative transfer calculations, and describe our so- 
lution technique. In our work we have assumed that the SN flow can 
be adequately described by a homologous expansion (i.e., v oc r). 
There are several reasons for doing this. First, SNe approach a ho- 
mologous state at late times. For example. Type la SNe rapidly ex- 
pand to a radius of < 10^'* cm from an object the size of the Earth 
(< 10^ cm) in one day. After such a rapid expansion, the assump- 
tion of a homologous expansion is excellent. This is also true for 
Type lb, Ic SNe which expand, from an object similar in size to 
the Sun, by over a factor of 100 in 1 day. For Type II SN the as- 
sumption becomes a good approximation after a few days for a 
blue-supergiant (BSG) progenitor and somewhat later (~ 20 d) for 
a red-supergiant (RSG) progenitor. Several factors contribute to the 
departure at early times. First the progenitor radius is large (par- 
ticularly for a RSG), and it takes time for the expanding ejecta to 
"forget" its initial radius. Second, the energy density of the radia- 
tion field at early times is significant compare d with the kinetic en- 
ergy, and this influences the hydrodynamics jFalk & Amett|[l977l ; 



^ In many papers it is referred to as the thermalization parameter, and can 
be thought of as the fraction of line photons that are "absorbed" on a given 
line interaction (e.g. [Hauschildt et al. 1992)- This probability not only al- 
lows for true absorption, but also, for example, re-emission of a photon 
from the upper level in another bound-bound transition. 



iDessart et al.ll2010bl) . Third, in Type II SNe, the reverse shock can 
set up a non-monotonic velocity field between the core and He/H 
interface. Finally, material can fall back towards the core, which is 
particularly important for high mass progenitors. 

The second reason for neglecting the departure from a homol- 
ogous expansion is that in this complex, but still exploratory work, 
we do not wish to be concerned with dynamics. A third reason is 
that the radiation transfer equation is simpler and there are compu- 
tational advantages in solving this equation compared with the full 
transfer equation. In particular, the third moment of the radiation 
field does not appear in the moment equations. 

This paper is organized as follows: In ij2]we briefly charac- 
terize how time dependence can affect the modeling of SNe. The 
radiative-transfer equation, its moments and associated boundary 
conditions in space and time, are discussed in SjS] The time depen- 
dent rate and radiative equilibrium equations are discussed in ij4] 
while the deposition of radioactive energy is discussed in !j5] A 
brief description of the numerical method used to handle the La- 
grangian derivative is provided in ij6l As in other CMFGEN model- 
ing, we use a partial linearization method to facilitate the simulta- 
neous solution of the radiative-transfer equation, and the rate and 
radiative equilibrium equations (f|7}. The use of super-levels, which 
reduces the size of the non-LTE model atoms, and problems intro- 
duced by the use of super-levels are discussed in ijs] Model con- 
vergence is discussed in i|9] while code testing is discussed in ijlOl 
The computation of the observed spectrum, which is done in both 
the comoving and observer's frames, is discussed in t]l II In m2\ 
we discuss the time dependent gray transfer equations, which can 
be used to provide an initial estimate of the temperature structure. 
From these equations we derive a global energy constraint (i |12.1| l. 
Finally, in i]13l we discuss additional work that is needed to im- 
prove ID modeling of SNe. 



2 TIME DEPENDENT EFFECTS IN SNe 

Time dependence enters into the radiative-transfer equation be- 
cause the speed of light is finite. To characterize its importance, 
we can consider three separate cases: 

(i) Consider an atmosphere with a scale height AR. In the opti- 
cally thin limit, the characteristic light travel time is simply AR/c. 
For a SN, the characteristic flow time is R/v and hence the ratio of 
the light travel time to the flow time is 

AR v 

which for v = 10, 000 km s"^ and AR/R = 1/10 (appropriate 
forp = po(''o/'")^") we have tiight/tflow ^ 1/300. Because this is 
much less than I, it is generally assumed that the time dependency 
of the radiative-transfer equation can be ignored when modeling 
the photosphere. This will generally be reasonable even when the 
scale height of the photosphere approaches the local radius. 

(ii) In the optically thick regime the effective light travel time 
will increase since photons cannot freely stream — instead they 
perform a spatial random walk, mediated by absorption and scatter- 
ing processes. In this random walk, photons will travel a distance 
AR from their origin on a timescale tAR/c, and hence 

_ tAR V 

flight /iflow — — ^ — — ■ 

That is, the random walk of the photons through the medium in- 
creases the light travel time by a factor of t. Thus, as soon as 
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the optical depths become large (e.g., r > 100) the random-walk 
time will approach, or exceed, the flow timescale. In this case, the 
time dependence terms are cmcial, and will control the tempera- 
ture structure of the SN ejecta. The preceding statement applies to 
all SN ejecta and their light curves, i.e., optical depth effects sys- 
tematically require a time-dependent treatment. However, the effect 
is most pronounced, and lasts longer, in those SNe (i.e.. Type IIP) 
with large ejecta masses. 

The random-walk time scale is equivalent to the diffusion time- 
scale, but the use of the word diffusion in this context is some- 
what of a misnomer, since "classical" diffusion noimally refers to 
the case where we have a temperature gradient. However, in core- 
collapse SNe the temperature gradient is imposed by the shock, and 
may be positive, negative, or zero. At large optical depths the pho- 
tons are trapped by the optically thick ejecta. If the SN were not 
expanding, some of these photons would eventually "diffuse" to 
the surface, and would impose a temperature gradient on the ejecta. 
However, in the "optically thick" regions of Type II SN ejecta the 
expansion time is shorter than the diffusion time, and the tempera- 
ture is set by the influence of expansion (adiabatic cooling) on the 
initial temperature structure, and photon escape from depth is pri- 
marily determined by the contraction of the photosphere to lower 
velocities. Initially the light curve reflects the cooling induced by 
adiabatic expansion; later the light curve is controlled by the es- 
cape of therma l energy as the recombination front moves into the 
SN ejecta (e.g.. lEastman et al.iri994 ; Amett 199 j) . 

The relative importance of the different energy sources will de- 
pend on the mass of ^^Ni created, and the mass and the initial radius 
of the ejecta. For example. Type la SNe differ significantly from 
Type II SNe. In Type la SNe the explosion occurs in such a small 
object (approximately the size of the Earth) that the initial explo- 
sion is never seerQ — rather the entire visible light curve of Type 
la SNe is dictated by the re-heating of the SN ejecta by '^'^Ni/^'^Co 
decay. In the unmixed Type Ilb/Ib/Ic models of lOessart et alj ^201 ih 
the re-brightening caused by the diffusion heat wave arising from 
the decay of ^^Ni was not seen until about 10 days — prior to this 
the light curve exhibited a plateau. Models without ^^Ni showed 
a similar plateau, but subsequently faded. The light curve of SN 
I987A, because it stems from the explosion of a BSG star, is ini- 
tially controlled by adiabatic cooling, but then, as in Type la SNe, 
beco mes dominated by the energy released frorn nucl ear decay 
(e.g.. lShigevama & Nomotolll99d ; lBlinnikov et alj|2000l) . 

(iii) A third case to consider is where the photospheric proper- 
ties are changing rapidly. Even if tught/iflow < 1, and light travel 
time effects are unimportant for the transport of radiation in the 
photosphere, light travel time effects can still influence observed 
spectra. This occurs because photons arising from an impact pa- 
rameter p = (i.e., the observer's line of sight striking the SN 
center of mass) will arrive at the observer at a time -Rphot /c ear- 
lier than photons arising from rays from an impact parameter with 
P ~ .Rphot (i-e., the observer's line of sight tangent to the photo- 
sphere). An example of this effect can occur near shock breakout 
where the photospheric conditions are changing rapidly. A RSG has 
a radius of ~ 10^'* cm, and hence the travel time across the radius of 



* iHoflich & SchaefeJ j20O9l) have studied whether the explosion could be 
seen in X-rays and gamma-rays. Their studies suggest that the Burst and 
Transient Source Experiment (BATSE) should have seen about 13 it 4 
nearby SNe la while none were detected. They suggest that absorption by 
the accretion disk can account for this discrepancy. The optical light curve, 
arising from the shock breakout, s hould be detectable, with an optical/U V 
luminosity of order lO*^ to lO'^ Lr^ \Y-ao et alj2010l ; lRabinak et al.ll201lh . 



the star is Ih which is comparable to the variability timescale, and 
hence light-trav el times effects w ould significantly effect the early 
fight curve te.g., lGezari et al.l2008 : Nakar & Sari 2010). While the 
above discussion is primarily related to early photospheric phases, 
it may also be relevant at late times. At late times, the evolution 
of the properties of the SN ejecta and light curve are controlled by 
the deposition of radioactive energy, and the relevant time scale to 
compare the light travel time with is the smaller of the decay time 
and the flow time. For the photospheric phase, the dominant nu- 
clear energy sources are the decays of *'^Ni and ^®Co, and these 
have decay times of the same order as the flow time. 

In the following, we will be primarily concerned with ef- 
fect (ii). Case (i) is difficult to handle because of the very short 
time scales. In our studies the effects are likely to be small, even 
allowing for the extension of the SN ejecta. While we do treat 
time dependence, the large time steps (2% to 10% of the current 
time) means that the (small) influence of freely traveling light is 
smoothed in time - in a 10% time step the light can travel the 
whole of the grid. To accurately treat such an effect, if important, 
we would need to consider time steps comparable to the light travel 
time across grid points. For our studies it is the flow time that is 
cmcial, since it is the primary factor determining the time scale on 
which the ejecta properties and radiation field change with time. In 
the limiting case of a very thin photosphere (i.e., the atmosphere 
can be treated in the plane-parallel approximation), case (iii) is 
more a bookkeeping exercise — the observed spectrum is simply 
the weighted sum of spectra from a series of steady-state plane- 
parallel atmospheres computed at different time steps. That is. 



2ty 

F^to+d/c) = — - 



pl„ {to - Az/C, Rm,vi,p)dp (I) 



with 



(2) 



where to is the time the observed light is emitted for an impact 
parameter p = J?max, d is the distance to the SN, and for simplicity 
we have ignored the expansion of the photosphere during the time 
interval At — Rma/c. 

In our simulations Rma is the same for all frequencies, and 
is chosen so that at the outer boundary the SN ejecta is optically 
thin. While the latter condition cannot always be met (e.g., in the 
continuum of the ground state of a dominant ion or for a strong res- 
onance transition) extensive tests with W-R and O star models, and 
a lesser number of tests with SN models, show that this does not 
affect the observed spectra in the important (i.e., the regions con- 
taining the flux) spectral regions (see i) |3.2| for further discussion). 
At early epochs, the outer boundary typically lies at a velocity of 
0.1 to 0.2 v/c. 



3 RADIATIVE TRANSFER 

The transfer equation in the comoving-frame, to first order in v/c, 
and assuming a homologous expansion (so that dv/dr = v/r) is 

1 dlu ^ + V dlv ^ (1 — /i^) dh 
c dt c dr r dji 

— 1 — " Xviv (3) 

rc ou rc 

jMihalas & Mihalad Il984h where = I{t,r, fi,^), = 
x{t, r, v) is the opacity, r\v — rj{t, r, v) is the emissivity, and /, x, 
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rj, jj,, and u are measured in the comoving-frame. rji, and Xi^ are as- 
sumed to be isotropic in the comoving frame, and we have dropped 
the subscript o (which denotes comoving frame) on I, x, V' ^ for 
simphcity. The moments of the radiation field are defined in terms 
of the specific intensity by 



[J, H, K] = J [1, ^, ^"]7(t, r, ^I, u)d^, . 



(4) 



Integrating over the transfer equation, and rearranging terms, we 
obtain the zeroth and first moment of the radiative-transfer equa- 
tion: 



1 D{r\J,) ^ 1 d{r^H,) 



Dt 



dr 



vv dJv 
rc du 



riu - XvJv (5) 



and 



-H — K 1 K — = ~Xvtlv 



Dt 



dr 



rc dv 



(6) 

where D/ Dt is the Lagrangian derivative which for spherical ge- 
ometry is 



D__ d_ d_ 
Dt ^ dt '^'"dr ' 
An alternative formulation of equations |5]and|6]is 

V duJv 



(7) 



1 D{r^J^) ^ 1 dir^'H^ 



Dt 



dr 



du 



= Vi' - XvJv (8) 



and 



1 D{r'^H^) 1 djr'^K^) - v duH^ 



-XvHjj . 
(9) 

These equations more directly show the connection to the gray 
problem (the frequency derivatives integrate to zero) (see ^\2\ and 
the importance of adiabatic cooling due to expansion in the opti- 
cally thick regime. That is J oc l/r'* and hence T cx 1/r. We 
utilized Eqs.[5}{6]in our formulation, primarily because, in the ab- 
sence of the D/Dt terms, they are identical to those we use to solve 
the transfer equation in stellar winds (at least if i; oc r). 



3.1 The Eddington factors, f„ 

The two moment equations discussed in ij3](Eqs.|5]&[6]l contain 3 
unknowns {J^, H,^, and A'^), and are thus not closed. To close these 
equations we i ntroduce the closure relation = K^/ Ji, (see, e.g., 
lMihalaslll978l) where is obtained from a formal solution of the 
transfer equation (Section [3. 3t . Since the formal solution depends 
on Sv, which in turn depends on J,j, we need to iterate. As is well 
known, convergence of this iteration procedure is rapid. 

The Eddington factor, which provides a measure of the 
isotropy of the radiation field is a ratio of moments. As such they 
are relatively insensitive to the radiative transfer, and we therefore 
compute the Eddington factors using the fully relativistic, but time 
independent, moment equations. By solving the time-dependent 
moment equations we coixectly allow for the trapping of radiation 
in the optically thick SN envelope. Thus the temperature structure 
is correct and we compute the correct level populations. Near the 
surface we are neglecting light-travel time effects, but this will only 
be important in the event that it effects the angular distribution of 
the radiation, and in particular the Eddington factor . 



3.2 Boundary Conditions 

3.2.1 Outer boundary 

For the solution of the moment equations we assume 

= KJ, (10) 

at the outer boundary where h is computed using the formal solu- 
tion. This assumption is used in Eq. [6] which is differenced in the 
usual way. The use of this boundary condition can sometimes gen- 
erate instabilities in the solution, particularly at longer wavelengths 
(in the infrared). One manifestation of this instability is that a nu- 
merical/differencing error introduced into the solution may remain 
constant (or even grow). Since J declines with increasing wave- 
length, the percentage en'or in J increases. The instability is mon- 
itored using the formal solution, and can be reduced/eliminated 
by choosing a finer resolution at the outer boundary such that 
Ari << Ar2 where indices, starting at 1 at Rmax, are incre- 
mented inwards until ND at the innermost grid point.. 

In the formal (i.e., ray) solution we have two choices at the 
outer boundary. In the simplest approach we adopt 1,7 = at 
the outer boundary. Depending on the size of the model grid, this 
choice can create model difficulties because > 1 at some fre- 
quencies implying that ~ Sv. Using = we obtain 
J,j — Sv = S,j/2 rather than Jv — Sv « Sv for these frequen- 
cies, and this causes rapid changes in level populations at the outer 
boundary. In order to treat the boundary region correctly it would 
be necessary to resolve the outer boundary using many additional 
grid points. In practice this is unwarranted — the outer boundary 
condition (provided it is reasonable) does not affect the solution, 
and the shaip truncation at the outer boundary is an artifact of the 
model construction. 

As an alternative, we can set I~ = at an extended outer 
boundary, and solve the transfer equation in this extended region to 
obtain at the true outer boundary. In SN models, the extended 
region typically extends a factor of 1.5 larger in radius, and opac- 
ities and emissivities in this region are obtained by extrapolation. 
In order to limit the velocity to be less th an c (the SN1987A mod- 
els discussed bv lDessart & Hillied ( [2010|) had u(i?max)/c = 0.205 
and (-Rext/c = 0.29) we use a beta-law extrapolation of v [v{r) = 
(1 — -Rcore/f)^] , rather than a pure homologous expansion. The 
technique of extrapolating the radius produces a much smoother 
behavior in the level populations at the outer boundary. It has been 
successfully used previ ously for the modeling of Wolf-Rayet stars , 
LBVs, and O stars (e.g.. lHillieJl987l : lHillier & Milleill998Lll999h . 

3.2.2 Inner boundary 

For core-collapse SNe, the inner boundary represents the innermost 
ejecta layer that sits at the junction between the ejecta and fallback 
material (deemed to fall onto the compact remnant). Depending on 
the model, and in particular the b inding energy of the progenitor he- 
lium core ( iDessart et al.ll20I0al) . the inner ejecta velocity can vary 
between <^ 100km s~^ for standard explosions while for energetic 
explosions, or explosions of low mass progenitors, the velocities 
can be an order of magnitude larger. For Type la SNe, for which 
no remnant is left behind, the inner core v elocities are on the order 
of > 1000 km s~^ jKhokhlov et alj 19931) , and still little mass trav- 
els at such low velocities compared to what obtains in Type II SN 
ejecta. 

For simplicity we adopt 



Hv = Q 



(11) 



© 201 1 RAS, MNRAS 000.[Tll20l 



6 



D. John Hillier, Luc Dessart 



(in actuality we tend to adopt a small non-zero value for H,j, but 
such that Hi, « J^). For models with an optically thick core, 
this assumption will produce a smooth run of level populations and 
temperature with depth near the inner boundary. Our assumption of 
H,^ = ignores the possible influence of magnetar radiation etc. In 
reality (and assuming we have the correct SN density structure) the 
flux at the inner boundary is determined via symmetry arguments. 
The inward intensity at the inner boundary is the outward intensity 
at the inner boundary redshifted in frequency space and from an 
earlier time step. Thus 

(12) 

with i/i = !/o(l + — and where we have ignored the 

distinction between ^ and jio since v « c. 

Thus Hv has a complicated frequency variation and would 
need to be determined in the iteration procedure. However, when 
integrated over frequency, H jSJ and thus from the gray per- 
spective the assumption of = is accurate. 

In practice it is reasonably easy to take the frequency shift into 
account. In the comoving-frame technique we iterate from blue to 
red, and thus the intensity information at the high frequency is al- 
ways available. In the formal solution we can either use //^ = 
(and thus = Iv) or we specify taking into account the fre- 
quency shift. Doing the later is of importance when the core is op- 
tically thin, and we are computing the observed spectrum. When 
computing an observed spectrum, and with an optically thin enve- 
lope, the frequency shift is potentially important for obtaining an 
accurate spectrum. 

We have also implemented the non-zero H,^ case into the 
moment solutions. The straight forward approach of using i/^ — 
hvJv at all frequencies introduced some numerical instabilities into 
the solution, and it was found that it was better to use this approach 
only in frequency regions where the core was optically thin. We 
also find some convergence difficulties at the inner boundary as 
we transition from optically thick to optically thin, and these are 
still being investigated. The cause of the difficulty is most likely 
related to resolution issues and the large variation in optical depth 
with frequency. When H,j is no longer small, we effectively have 
a "photosphere" at the inner boundary, and in this "photospheric" 
region populations will change rapidly. 

Accounting for the time difference in the inner boundary con- 
dition is more problematic (since the intensities at earlier time steps 
would be needed) but its effect is likely to be small since the rele- 
vant time scale is very small (At < 2i?„iin/c = 2tage {v{Rmin) / c) 
where iage is the time since the SN exploded). Theoretically, it 
could be taken into account using I_nu computed at previous time 
steps. 

Our adopted boundary condition is physically realistic, and is 
a much better approach than applying a boundary condition at some 
intermediate depth below the photosphere. It also allows a better 
handling of the spectral evolution at late times when some spectral 
regions start becoming transparent while others remain thick. 

3.2.3 Time 

The radiative transfer problem represented by Eqs. |3] |5] & |6]is 
an initial value problem in time. Thus, to solve these equations, 
we need to fully specify {tinu , r, /i) and the coiTesponding mo- 
ments. The computation of Iv{tinit,r, /i) can only be done in con- 



junction with a full hydrodynamical simulation. This is beyond the 
scope of the present work. We therefore proceed as follows: 

(i) We adopt the temperature, density, and abundance structure 
from a full hydrodynamical model. In the past, Stan Woosley has 
supplied us with such models. As these models, which are for RSG 
progenitors, are not fully resolved in the outer layers, we stitch a 
power law onto the hydrodynamical models at some suitably cho- 
sen velocity (e.g., 10,000 to 30,000 kms~^). As our choice may 
affect spectra, particularly at early times, we now obtain/compute 
models that more accurately describe the outer layers. The issue is 
less of a problem for Type la, lb, and Ic SNe after one day, because 
the outermost layers of these models have been greatly diluted by 
the huge expansion the ejecta has undergone. 

(ii) We solve for the radiation fields using a fully relativistic 
transport solver ignoring time dependent terms (dJ^/dt, dH^/dt). 
In performing this calculation we allow for non-LTE, and we iter- 
ate to obtain the populations simultaneously with the radiation field 
(i.e., we undertake a normal non-time dependent CMFGEN calcula- 
tion). The temperature is held fixed at the value obtained from the 
hydrodynamical model. 

At depth this procedure is more than adequate for Type II SN, since 
,J,j is close to Bv. However closer to the photosphere, effects may 
be seen since departs from B^, and the diffusion time can be 
significant relative to the age of the SN. The long term effects of 
this inconsistency are probably small but need to be examined. We 
can examine errors introduced by the starting solution by compar- 
ing CMFGEN results at some future time with the hydrodynamic 
solution at that time, and by comparing models at the same epoch 
but which begin with different starting solutions. 

A second issue is that the temperature structure from the hy- 
drodynamical structure will not be exact — typically the temper- 
ature has been determined assuming LTE and a gray-opacity ap- 
proximation. Only on the second time step will the temperature 
structure of the outer layers adjust itself to be more fully consis- 
tent with our more accurate treatment of the radiation field, and the 
non-LTE opacities and emissivities. Since the LTE approximation 
is valid for the inner layers, inconsistencies in the outer tempera- 
ture structure will only affect the earlier model spectra — with time 
these outer layers become optically thin and no longer contribute to 
the observed spectrum. A comparison of the different temperature 
structures f or Type Ib/Ic SN com puted using different assumptions 
is shown bv lDessart et al.l ( l20Ill) . 

For the time step, we typically adopt At = O.lt. Tests show 
this is adequate over most of the time evolution. Exceptions will 
occur at shock breakout (not currently modeled) and at the end of 
the plateau phase in Type II SN as the photosphere passes from the 
hydrogen rich envelope into the He core. As the plateau stage ends, 
the luminosity can change by a factor of a few on a time scale of 
10 days. As this is roughly O.I times the age of the SN, a smaller 
st ep size needs to be used. In the Type II-P core collapse models 
of lOessart & Hillieil ( |20I ih time steps of 2 to 5 days did a much 
better job of resolving t he end of the plateau phase (see Fig. 10 in 
iDessart & Hillieill201 ih . At present the time step is set by hand, 
with 10% being the default. 

3.3 Formal Solution 

In order to solve for the moments of the radiation field we need 
to compute the Eddington factors, which are used to close the mo- 
ment equations. These are computed using the fully relativistic, but 
time independent, transfer equation. The fully relativistic transfer 
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equation is 

7(1 + M di 



dl 



+ 7(1 -m) 



-A 



21 

dfj. 



+ 37 



/3(1-m') 



/3(1-m') 



+ AtA 



+ M 



a/ 

I = V-Xl (13) 



where jS = v/c. 



7 = l/\/l-/32 



(14) 



(15) 



and it is explicitly assumed that J = I(t, r, ^, t /) with /, fi and 
measured in the comoving frame jMihalasll 19801) . 

To solve the transfer equation we will assume that we can 
ignore all time-dependent terms. This is reasonable, since we are 
not using this equation to compute the temperature structure of the 
ejecta. Rather, we will use its solution to compute the Eddington 
factors for the solution of the moment equations. As discussed in 
!j2]and ^3.U errors introduced by this approximation are likely to 
be small. For the computation of the observed spectrum, we are ig- 
noring the light travel time for freely flowing photons, but the errors 
introduced by this assumption will generally be small, except near 
shock breakout (Secti on [2). The techn ique we use to solve Eq. [T3] 
closely follows that of lMihalasI ( Il980l) . 

Using the method of characteristics, the transfer equation can 
be reduced from a partial differential equation in three dependent 
variables (r, v, fx) to a series of equations in two dependent vari- 
ables (s, v). From Eq.[T3]it follows that the characteristic equations 



the current frequency, we have 



ds 



= 7?fc + a-klk-i - ixk + 311-1- ak) Ik 



with 



t^k 



-n 



I'fe-i - i^k 

This may alternatively be written as 

dh 



dTk 



= Ik ~ Sk 



with 



and 



dfk 



■ (xfe + 3n + ak) ds 



Sk = iVk + akh-i)/ iXk + 3n + Ofe) 



(21) 



(22) 



(23) 



(24) 



(25) 



Since Eq.|23]is the standard form of the transfer equation along a 
ray it can be solved by the usual techniques. Our approach is to 
assume that Sk can be represented by a parabola in Tk over three 
successive grid points, and we analytically integrate the transfer 
equation using this representation. We first integrate inwards to ob- 
tain I^ , then outwards to obtain I^ . The radiation moments are 
obtained by numerical quadrature assuming is a linear function 
of ^ — this is preferred over higher accuracy schemes because of 
stability. 

More recently an alternative approach has been suggested. In 
this mixed frame approach the transfer equation is written in terms 
of an affine parameter, and the como ving frequency (or wavelength; 
IChen et al.l2007l : lBaron et al.l20d^ . Its advantage is that the trans- 
fer equation is solved along geodesies (straight lines in our models) 
allowing the approa ch to be easily g eneralized to three dimensions. 
Like the method of lMihalasI ( ll98ol) . the opacities and emissivities 
are evaluated in the comoving-frame. 



dr 
ds 



= 7(a^ + P) 



and 



dH , 2^ 



i + Z^M 



A 



(16) 



(17) 



We choose Nc rays to strike the stellar core, while the remain- 
ing rays are set by the impact parameter "p" which is defined by 
the radius grid. To obtain the curved characteristics we integrate 
the above equations using a Runge-Kutta technique. The accuracy 
of the integration can be checked by using the relation between /is 
and /io, and noting that [is is simply \/[l — (p/r)^]. 

As an aside we note that the dp/dt term may be written as 
Dp/Dt — vdP/dr. Thus for a Hubble flow A reduces to 



A ^ fi 



dl3_ 
dr 



(18) 



Along a characteristic ray the transfer equation reduces to 
dl dl 

— ~i.n— =7?-(x + 3n)/ (19) 
OS av 



with 



n = 7 



/5(l-A'') 



/iA 



(20) 



Using backward differencing in frequency, and using k to denote 



4 TIME DEPENDENT RATE AND RADIATIVE 
EQUILIBRIUM EQUATIONS 

The time-dependent rate equations can be written in the form 



Dn,/p _ 1 D{r^ni) 



Dt 



Dt 



TliRi 



(26) 



where p is the mass density, is the population density of state i 
for some species, and Rji is the rate from state j to state i. In gen- 
eral, the Rji are functions of Ue (the electron density), temperature, 
and the radiation field (see, e.g., Mihalas 1978), and when charge 
exchange reactions are included they are also dependent on Ui. 
Thus the equations are, in general, explicitly non-linear, although 
with the exception of the temperature, this non-linearity is only 
quadratic. Of course, there is a strong implicit non-linearity result- 
ing from the dependence of the radiation field on the populations. 
In CMFGEN we allow for the following processes: bound-bound 
processes (including two photon decay in H and He I), bound-free 
processes, collisional ionization and recombination, collisional ex- 
citation and de-excitation, charge exchange (primarily with H I and 
He I), and Auger ionization ( Hillier. 1987;.Hillier & Miller.1998.) . 
The energy equation has the form 



De P Dp r ^ , ^ , . 



(27) 
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where e is the internal energy per unit mass, P is the gas pres- 
sure, and Edccay is the energy absorbed per second per unit vol- 
ume. When gamma-ray deposition is assumed to be local (see Ap- 
pendix), it is simply the energy emitted by the decay of unsta- 
ble nuclei (excluding the energy carried off by neutrinos). To im- 
prove numerical conditioning, continuum scattering terms are ex- 
plicitly omitted from the absorption and emission coefficient in 
equation|27] e can be written in the form 



where 



and 



eK 



eK + ei , 



3kT{n + ne) 
2jj,m,n 



ei 



E 



UiEi 
fimn 



(28) 



(29) 



(30) 



In the above, n is the total particle density (excluding electrons), 
m is the atomic mass unit, fi the mean atomic weight, T is the gas 
temperature, and Ei is the total energy (excitation and ionization) 
of state i. 

The influence of the time dependent terms in the rate and 
radiative equilibrium equations on eject a properties and syntheti c 
spectra has previously been discussed by [Dessart &Hillieil ( 120081) . 
They found that these terms are crucial for explaining the strength 
of the Hi lines in Type II-P spectra during the recombination 
epoch. Importantly, the terms also increase the He I line strengths 
in early spectra so that their line strengths are in better agree- 
ment with observation. However the terms do not just influence 
H and He lines — other lines are also affected since the elec- 



tron density chang es l lDessart&Hillieill2008l) . As noted earlier, 
lUtrobin & Chugaf feOOSi) . showed that the time dependent terms 
were also important for modeling the H I lines in the Type II- 
peculiar SN 1987A. More detailed work on the influence of 
time-dependence on Ty pe IIP SN spectra has been presented by 
lDessart & HillieJfeOlOl) . 



5 DEPOSITION OF RADIOACTIVE ENERGY 

During a SN explosion radioactive elements are created (for a 
review, see, e.g., IXrnett 199^. Type II SNe produce between 
0.01 and 0.3 M© of ^'^Ni (e.g.. iHamuvl |2003|) while a Type la 
SN can eject anywhere between 0.1 to 1.0 Mq of ^^Ni (e.g., 
iHoeflich & Khokhlo^^[l993 : IStritzinger et alj|200d) . The influence 
of the radioactive material on the light curve (and spectrum) very 
much depends on the ^''Ni mass and distribution, the size of the 
progenitor, and the ejecta mass. Radioactive decay powers the light 
curve of all SN types during the nebular phase (excluding those 
powered by an external source such as a magnetar). The decay of 
^^Ni (and its daughter nucleus, ''^Co) powers the entire light curve 
of Type la SNe. For SN 1987A, associated with the explosion of 
a BSG progenitor star, ^®Ni heating begins to become apparent at 
approximately day 10, and is responsible f or the broad maximum, 
which lasts for approxima tely 100 days dShigevama & Nomotol 
ll99(]| : lBlinnikov et alj|2000l) . 

The most important nuclear decays are 



^'^Ni '"^Co + j (half life = 6.075 days, 

AE = 1.7183 MeV per decay) 

and 

^'^Co "''Fe+7 + e+(half life = 77.23 days, 

Ai5(7-rays) — 3.633 MeV per decay; 
AiJ(positrons) — 0.1159 MeV per decay). 

The energies and lifetimes are all from 
[http://www.nndc.bnl. gOv/chart/l which for ^''Ni and ^®Co are 
based on the work of fauo et al.l l ll987l) . Energy release in the form 
of neutrinos is neglected, since ejecta densities at the times of 
interest are all well below the neutrino-trapping densities. 
To implement the decay of the radioactive elements we analytically 
solved the 2-step decay chains. As a consequence, the accuracy of 
time-dependent abundances is not affected by our choice of time- 
step. Another corollary of our approach is that we do not use the 
instantaneous rate of energy deposition - rather we use the aver- 
age rate of deposition over the time step. In practice, the difference 
between the two values is small. Our nuclear-decay reactions are 
currently being updated so that we can accurately follow the evolu- 
tion of the SN ejecta abundances. 

In Type II SNe it is reasonable to assume that all the energy, 
emitted as 7-rays or positrons, is deposited locally. At later times 
(but well after maximum brightness) this assumption breaks down, 
and 7-rays can travel long distances depositing energy elsewhere 
in the SN, or escaping entirely. On the other hand, in Type I SN, 
the small ejecta masses mean that non-local 7-ray energy deposi- 
tion, and 7-ray escape, especially in models with extensive mixing, 
needs to be taken i nto account even at relatively early times (e.g., 
iHoflich etalj 19921) . We have developed a 7-ray transport code (see 
Appendix), to allow for non-local 7-ray energy deposition. 

At present we simply assume that the energy from radioactive 
decays goes into the thermal pool. This is likely to be valid for mod- 
eling the early light curve of most Type II SNe (but well beyond 
maximum brightness) but may fail relatively early for other SN^fl 
One factor limiting the influenc e of non-thermal proc esses at earlier 
times is the mass of the ejecta toessart et al1l2012bl) . Type IIP SN 
have much greater ejecta masses than do most Type I SN, and hence 
the photosphere at early times is confined to the outer layers where 
^®Ni has not been mixed. Non- thermal ionization and excitation 
processes, for example, are believ ed to be important for explainin g 
the appearance of He I in SNe lb jLucvlll99ll ; lDessart et al.ll201 ll) . 
and for the production o f H I lines in the nebular phase of Type II 
SNe and in SN 1987A l IXu & McCra"vlll99ll : lKozma & Franssor] 
Il992lll998allbl) . Non-thermal excitation and no n-thermal proce sses 
have recently been incorporated into CMFGEN jLi et alj|2012l) and 
we are currently inves tigating their infl uence on the spectra o f SN 



we are currently mves tigatmg their intl uence on tne spectra o r aJN 
I987A jLi et alj|2012h. Type Ibc SNe foessart et alj|2012bl) , and 
Type la SNe jDessart et al.ll2012al) . 



^ Type II SNe have much larger ejecta masses than do Type I SNe. As a 
consequence radioactive material is generally hidden from view for much 
longer. Deep below the photosphere, conditions ai'e highly ionized and close 
to LTE, and any radioactive energy deposited will be thennalized. 
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6 NUMERICAL SOLUTION 

To handle the Lagrangian derivative we adopt an implicit approach 
( iHoflich et all 1 19931 : fstone et all 1 19921) which is explicitly stable, 
for both the radiative-transfer equation, the rate equations, and the 
radiative energy balance equation. Thus we write 



fDX_\ _ Xi- 



(31) 



\ dt J ^ ti — ti-i 

where i refers to the current time step, and i~ 1 refers to the previ- 
ous time step. As our approach is only first order accurate in time, 
the time step must be kept "small". An explicit approach cannot 
be used as the conditions on the time step would be unnecessarily 
restrictive (e.g., IStoneetal.lll992h . 

All other terms in the radiative-transfer equation and rate 
equations are evaluated at the current time step. Since the quantities 
at the previous time step are known, the Lagrangian derivative sim- 
ply introduces extra source terms into the equations. With the above 
scheme the zero-order moment equation, for example, becomes 

1 dr^H^\ fuvdJ^' 
rc dv ^ 

3 

Mi + 



dr 



Ti 



cAt 



(32) 



This has exactly the same form as the equation nor mally solved in 
CMFGEN, a nd can be differenced in the usual way jMi halas et al.l 
Il975lll976l) . The only difference is that there is a modification to 
the opacity and emissivity arising from the radiation field at the pre- 
vious frequency. The solution of the rate e quations has previously 
been discussed bv lDessart & Hilliej Jioosl) . 

An alternative approach would be to choose a semi-implicit 
approach with all terms evaluated at the midpoints of the time step. 
This would be second order accurate, and should also be stable. 
Its disadvantage is an increase in complexity — ad ditional terms 
would need to be saved from the previous time step. iHoflich et al.l 
( Il993h use an adjustable differencing approach, controlled by a pa- 
rameter 9, which allows for an implicit approach, a semi-implicit 
approach, or a fully explicit approac h. Another p ossibility, utilized 
in some evolutionary codes and by ILucvI ( |2005|) , is to generate a 
difference formula using the two previous time steps. 

At depth the opacity and emissivity are large, and generally 
the time terms in the "modified opacity" and "modified emissivity" 
are almost negligible (unless a very small time step is used), and 
their influence on the radiation transfer at a given frequency could 
be neglected if one were only interested in the formal solution of 
the transfer equation. However the terms are of crucial importance 
for the energy balance, since their importance is enhanced when we 
integrate the transfer equation over frequency ( iil2t . In particular, 
if we assume radiative equilibrium holds the opacity and emissivity 
are "analytically" cancelled from the integrated transfer equation, 
and thus the time dependent terms can be very important. As can 
be seen from the gray transfer equations (i ]12t . the temperature at 
large optical depths is almost set entirely by the time-dependent 
term, and for a homologous flow without energy deposition from 
radioactive decay it simply scales as 1/t. 

For most of our calculations we have set the time step to be 
10% of the current SN age (Section 13. 2. 3) . At some epochs, a 
smaller time step may be needed. Since we need to iteratively solve 
the radiative-transfer equations in conjunction with the rate equa- 
tions, a factor of 2 reduction in the time step does not mean a factor 



of 2 extra computational effort, since each time step will require 
fewer iterations (because the starting solution will be closer to the 
final solution). 

As we are considering SN ejecta in homologous expansion, 
we use the velocity as our Lagrangian coordinate. 

Before beginning our simulations, the hydrodynamical mod- 
els are mapped onto a revised grid optimized for the computation of 
the radiation transfer. The new grid is defined in terms of Rosseland 
optical depth, and generally we require at least 5 points per decade 
of optical depth. Similar grids are generally used in the modeling of 
stellar atmospheres, and are preferred to a grid defined in mass and 
optimized for the hydrodynamic calculations. At present we gener- 
ally perform a new mapping at each time step (although the map- 
ping is done using the calculations at the previous time step rather 
than a hydrodynamical model). Two problems can potentially arise 
with our mapping scheme: 

(i) In standard ID SN models there are often rapid changes in 
spatial composition and density distribution. Due to the numerous 
re-mappings, these rapid changes suffer from numerical diffusion 
in space. At present we do not make any special attempt to limit 
this diffusion. Numerical diffusion is less of a problem in mixed 
models because the spatial gradients are smaller. The steep density 
gradients are also stronger in core-collapse SN ejecta because of 
the chemically-stratified progenitor envelope structure. 

(ii) In some models (particularly for SN 1987 A) ionization 
fronts can arise. Because the optical depth (e.g., in the Lyman con- 
tinuum) can vary rapidly across the front, complications arise in 
the radiative transfer and in model convergence. To overcome these 
problems we semi-automatically revise the grid across the front. 
Crudely, we adjust the grid so the front is better resolved. A con- 
trol file, which can be altered while CMFGEN is running, is used to 
specify parameters for adjusting the grid. If necessary, the grid can 
be adjusted at each iteration. 



7 LINEARIZATION 

The transfer, energy, and rate equations are a coupled set of non- 
linear equations. Thus their simult aneous solutio n requires an it- 
erative technique. As described by lHillieij ( ll990l) we adopt a lin- 
earization approach. While cumbersome to implement, we have 
found this technique works extremely well for the modeling of 
hot stars and their stellar winds. Another common approach is to 
utilize approximate lambda operators to allow for t he coupling of 
the lev el populations with the radiation field (e.g., IWemeil [19861 : 
iLanz & Hubenv 1995; Hauschildt & Baron 19 9j) . In so me cases 
this is combined with other approaches — iHoflichl ( l2003h also uti- 
lizes an equivalent-two level-atom approach in his SN modeling. 

The linearization of the radiative-transfer equation and the rate 
equations is straightforward, and as it is very similar to that pre- 
sented bytHillier (199(3) it will be described only briefly. In our 
approach we linearize both the rate equations and the radiative- 
transfer equation. Using the linearized radiative-transfer equation 
we can solve for the 5Jy in terms of 5xi^ and Srj,^, and hence in 
terms of Sni, 5n2, Sn-j,, . . . SNe, and ST. In principle, the 5J^ de- 
pend on the perturbed populations at every depth, however we only 
allow for the coupling locally (diagonal operator) or between adja- 
cent depths (tridiagonal operator). Using these equations we elimi- 
nate the radiation field (SJ,^) from the linearized rate equations. For 
the diagonal operator we obtain a set of Nv simultaneous equations 
(where Nv is the total number of variables) at each depth. For the 
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tridiagonal operator we obtain a block tridiagonal system of equa- 
tions (Nd X Nv in total, where Nd is the number of grid points). The 
latter can be efficiently solved using a block-matrix implementation 
of the Thomas algorithm for Gaussian elimination. In a typical SN 
model Nd = 100 and Nv ~ 2000, yielding 200 000 equations. The 
2000 super-levels typically represent 5000 to 10000 atomic levels. 

Starting estimates for the temperature and level populations 
are obtained as follows. For the temperature structure we either 
adopt a scaled version of the gray temperature structure (see i|12b . 
or simply adopt the temperature from the previous solution. For the 
level populations we adopt the departure coefficient from the previ- 
ous time step, although we also have the option to adopt LTE values 
(useful only for the first model of a time sequence). 

To solve the linearize rate equations we use LU decomposition 
followed by backward substitution utilizing standard LAPACK0 
routines. As discussed by .Hilliei. C2003,) . we improve stability by 
doing the following: 

(i) We scale the simultaneous equations so that we solve for the 
fractional correction Srii/ni. 

(ii) We use the LAPACK routine DGEEQU to pre-condition the 
matrix. 

These two procedures work very well, and we (generally) have no 
difficulty in solving systems of equations with Nv= 2000. When 
problems do occur it is usually because of poor starting estimates. 

The advantage of the diagonal operator is that the construc- 
tion of the linearization matrix is faster; each depth can be treated 
independently; in the early stages of the iteration procedure con- 
vergence can be more stable; and it requires less memory (nomi- 
nally down by a factor of 3 but because of our implementation we 
save only a factor of 2). The advantage of the tridiagonal opera- 
tor is its speed of convergence — far fewer iterations are needed 
to reach a specific convergence criterion. With both operators, Ng 
acceleration (Ng 1974; Auer 1987) can be used to accelerate the 
convergence. 

The time taken per iteration is very model dependent, and 
depends on the type of iteration. Iteration steps (for Nd=125, 
Nv=2319, ~400,000 bound-bound transitions) in which the 
linearization matrix is held fixed take about ~15 minutes, 
A— iterations take somewhat longer (~20 minutes) while a full it- 
eration may take 4 hours. The computations were formed using 
4 processors (Quad-Core AMD Opteron[tm] Processor 8378; 800 
MHz), and the computational times are about a factor of 2 longer 
than that obtained using a Mac Pro. The computational time scales 
as Nd^ 



8 SUPER-LEVELS 

To reduce the complexity of the non-LTE problem it is usual to 
use super-levels (SLs), or a variant thereof, to reduce the num- 
ber of variables, and hence the number of rate equations that 
need to be solved (e.g., Anderson 'l989';'Dreiz ler & Wemeij|l993l ; 
iHubenv & Lanz 1995; Hillier & Miller 1998). In the procedure 
adopted by iHillier & Milled ( |1998|) . SLs are used as a means of 
reducing the number of level populations that need to be deter- 
mined, while all transitions are still treated at their correct wave- 
length. Within a SL we assume that the states have the same depar- 
ture from LTE (although more complicated assumptions can also 

^ LAPACK is a software library for numerical linear algebra. Documenta- 
tion and the routines are available at|http://www.netlib.org/lapack/, 



be adopted); thus this procedure can be considered to be the natural 
extension of LTE. 

However SLs can introduce inconsistencies which may ad- 
versely affect convergence. Consider a three level atom, with two 
excited levels which are relatively close in energy, and which are 
treated as a single SL in the rate equations. 

The rate equation for the SL contains terms of the form 

712^21-^21 + 713^31^31 

while the radiative equilibrium equation will contain terms of the 
form 

7/21712^21.^21 + ft 7^31713^31^31 

where Zij is the net radiative bracket, Aij is the Einstein A coef- 
ficient, h is Planck's constant, and Vij is the transition frequency. 
When levels 2 and 3 are part of the same SL, and if the bound- 
bound transitions are the dominant transitions determining the pop- 
ulation of the SL and are a significant coolant, an inconsistency 
arises since the ratio 712/113 is fixed by the assumption that two 
levels have the same departure from LTE. In extreme cases this in- 
consistency can ca use models not to converge. 

As noted by iHillieJ | |2003|) , another manifestation of this 
inconsistency is that when radiative equilibrium is obtained the 
electron-energy balance equation will not be satisfied. The level 
of inconsistency gives an indication of the e rrors introduc ed by 
the use of SLs. To overcome the inconsistencv, iHillieJ ( l2003b sim- 
ply replaced 7/31 and U21 by an average frequency, v in the radia- 
tive equilibrium equation. This procedure overcame convergence 
difficulties and yielded better consistency of the electron-energy 
balance equation with the radiative equilibrium equation. Detailed 
tests, performed by splitting the important SLs, confirmed the va- 
lidity of the procedure in the atmosphere/winds of massive stars. At 
high densities the scaling approximation can be switched off. 

In SNe the simple scaling described above has potential dif- 
ficulties. Firstly, the lower densities and higher velocities increase 
the importance of line scattering compared with continuous pro- 
cesses. Second, the radiative equilibrium equation appears on the 
right-hand side of the zeroth moment of the gray equation. Since 
the above scaling is not done in the transfer equation there is an 
inconsistency between the transfer equation and the energy balance 
equation. In the present case the assumption of radiative equilib- 
rium is untrue, due to the work done on the gas, and the deposition 
of energy from nuclear decays, but there will still be an inconsis- 
tency between the transfer equation and energy balance equation. 

To remove the inconsistency, but at the same time retain the 
benefits of the SL approach, we now simply scale the line opacity 
and emissivity. The emissivity and opacity due to the transition I- 
3 is scaled by u/ph - it is necessary to scale both the emissivity 
and opacity so that B — i;/x at depth. Most of the scalings are 
less than 10%, and thus of the same order of accuracy as metal line 
transition probabilities. In principle the influence of the assumption 
is easily tested by altering the SL assignments - in practice this can 
be difficult because of memory requirements. Alternatively, after 
the temperature has converged, the scaling can be switched off, the 
temperature can be held fixed, and the populations can be updated. 
With this technique the direct influence of the scaling on the ob- 
served spectrum can be removed. No scaling is used when comput- 
ing observed spectra with CMF_FLUX ( ^TTFl 

^ As noted by the referee it is possible to define an average A and Z value 
for a transition between two super-levels. In practice this would remove the 
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9 CONVERGENCE 

In type II SNe we expect the temperature at large optical depths 
to evolve adiabaticalljij as l/r from the initial configuration, and 
as this is a local scaling we expected convergence to be achieved 
despite the extremely high optical depths. However, initial mod- 
els failed to converge, with the temperature stabilizing at the gray 
temperature adopted for the initial estimate. Convergence was only 
an issue with the temperature structure. When the temperature was 
held fixed, populations would converge rapidly at all depths. 

The lack of convergence was traced to two issues, and both 
are inherently linked to the large optical depths. Because of the 
large optical depths, any term included in the linearization of the 
energy balance equation must be accurately linearized in order that 
the correct cancellation can occur. 

Consider the following contribution, by a single line at a single 
frequency, to the radiative energy balance equation. 



-3 



= {xc + (j)xi) J ~ {lie + (t)-ni) 



(33) 
(34) 



In the above, the subscripts "c" & 'T' are used to denote the con- 
tinuous and line contributions, and <j) is the intrinsic line emis- 
sion/absorption profile. As the medium is optically thick we have 
J = S and hence 



J 



+ e<t)Si 



1 + 1 



(35) 



(36) 



with 6 = Xi/Xc- Obviously x{J — S) = 0. However, if Sc ^ 
5"; the line and continuum term will make contributions that are 
identical in magnitude, but of opposite sign. This cancellation must 
also be maintained in the linearization. The continuum contribution 
XciJ — Sc) to xiJ — S) is, for example. 



In our original version of CMFGEN we integrated the line con- 
tribution to the linearization over the full line before adding it to 
the linearized energy balance equation. For simplicity, and compu- 
tational requirements, we added the linearized continuum contribu- 
tion to this term using data at the final frequency in the line. How- 
ever, as the continuum varies across the line, full cancellation was 
not achieved. Since cancellation was not achieved we effectively 
overestimated the effect of any change in the temperature on the 
energy balance. This, in turn, meant that the corrections needed to 
achieve balance were severely underestimated, and hence stabiliza- 
tion of the temperature structure occurred. The problem was solved 
by performing the continuum linearization at each frequency in the 
line (although this is only necessary, and is only done, for the en- 
ergy balance equation). 

A similar problem arose from the use of impurity levels. The 
concept of impurity levels was introduced into CMFGEN to mini- 
mize the number of levels that were fully linearized. The levels are 



inconsistency when the radiative equilibrium equation is defined in terms 
of the net rates although in practice it will yield results similar to our fre- 
quency scaling procedure. However this procedure does not solve the prob- 
lem with the radiative equilibrium equation on the RHS of the transfer equa- 
tion where there is no distinction between Unes and continuum. 
* Radioactive decay will modify this simple scaling. 
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Figure 1. Illustration of the corrections to the temperature at a depth 
Tioss = 1 ■ 5 X 10*' in a Type II SN model. The model, at t ime 2.3 d after the 
explosion, is similar to those used to model SN 1987A bv lDessart & Hillieil 
(2008, 2010). The convergence shown is fairly typical of the convergence 
behavior at depth. The temperature correction for iterations 1 to 39 is zero. 
All models begin with A iterations until the corrections stabilize, and then 
a mixture of lambda and non-A iterations with T fixed at all depths is per- 
formed. An Ng acceleration was performed after iteration 93 (hence the 
spike at iteration 94). 



linearized for rates associated with their own species, but their in- 
fluence on other species is neglected in the linearization. The use of 
impurity levels does not affect the solution; it only has an influence 
on the convergence of the code. In our implementation, exact can- 
cellation was not achieved for terms involving impurity levels, and 
the rate of convergence was severely curtailed in the inner region. 

To illustrate the convergent properties of the models we illus- 
trate in Figs.[T}j2]the convergence properties of a typical Type II SN 
model. Only the temperature is shown. In general, it is the tempera- 
ture that determines the global convergence properties of the model 
— when the temperature is held fixed convergence is usually rapid. 
For the model shown convergence is achieved in 100 full iterations, 
with an iteration defined as any step which changes the populations 
(e.g., a A iteration, a full iteration, or an Ng acceleration). 

While the temperature corrections in Fig. [T] are small, they 
help to illustrate that our linearization procedure does achieve con- 
vergence at large optical depth. As noted above, in early work tem- 
perature corrections stabilised, and it was unclear whether conver- 
gence was being achieved. When we use the gray temperature so- 
lution at depth for the initial temperature estimate, the predicted 
temperature corrections are small. However, when we use the tem- 
perature from the previous time step, the error in T is of order 10% 
(for a 10% time step) and we find good convergence provided we 
keep the populations consistent with the current temperature struc- 
ture. It also should be noted that the actual error achieved at a given 
iteration depends on both the size of the correction, and the rate at 
which the correction is changing. In Fig. 1, the true error, as es- 
timated over iterations 55 to 60, is approximately a factor of 10 
larger than the current correction. In Type II SN it is the shaip H 
ionization front, which moves to smaller velocities as the SN ages, 
that controls convergence of the model. 

Corrections to the populations are generally much larger than 
the temperature coiTections, and tend to behave more erratically, es- 
pecially in the early iterations. Experience has allowed us to anive 
at the following iterative procedure: 
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Figure 2. Same as Fig.[T] but now showing the convergence at a depth tj^q^s of ~ 1000 (left), ~ 1 (middle), and ~ 0.001 (right). The two types of symbols 
are used to indicate the sign of the correction. 



(i) Perform A iterations until the maximum correction is less 
than 50%. The temperature is held fixed. 

(ii) Perform full iterations, with the temperature held fixed, until 
the maximum correction is less than 50%. If the maximum correc- 
tion is greater than 200%, we revert to A iterations. 

(iii) When the maximum correction on a full iteration is again 
less than 50%, allow the temperature to vary. If the maximum cor- 
rection is greater than 200%, we revert to A iterations. When we are 
no longer switching between full and A iterations, we only compute 
the variation matrix (i.e., the linearized rate equations) every third 
iteration. 

(iv) When the maximum correction is less than ~ x% {x is an 
adjustable parameter typically set to 5) we hold the variation matrix 
fixed. It is recomputed if the maximum correction becomes larger 
than 3 x x%. 

(v) When 4 successive iterations are less than VAL_DO_NG 
(typically a few percent) we perform an Ng acceleration. We re- 
peat Ng accelerations every 10 to 20 iterations. 

This procedure is automatically handled by CMFGEN, al- 
though it is also possible to manually force a different behavior by 
restarting the code. Many of the control parameters are adjustable, 
although in practice we tend not to change them very much. 

In hot-star models we also have the option to fix the temper- 
ature automatically in the outer region until the populations have 
stabilized. This generally works well in hot-star winds, but in SN 
models it caused convergence problems around ionization fronts. 

A typical model is stopped when the maximum convergence 
is less than 0.1%, although most levels have already achieved a 
higher convergence. Several output files can be checked to verify 
that true convergence has been achieved. Some models in a se- 
quence are converged to a much higher accuracy (e.g., 0.001%) to 
further prove convergence. Convergence to higher accuracy is gen- 
erally not time consuming since the variation matrix can generally 
be held fixed — it is the initial iterations, with a mixture of A and 
full iterations, that take the most time. 

In Type II SN models two factors can affect the speed of con- 
vergence. 

(1) At a few depths (often only a single depth), a few lev- 
els may oscillate and fail to converge. These levels are generally 
unimportant, and do not affect the spectrum. However, they can 
affect convergence since it is the maximum change that controls 
the iteration procedure and ultimately determines the convergence. 
NG accelerations, averaging the last two full iterations, and forcing 
multiple A-iterations, can overcome the problems. More recently 
we have implemented a procedure to use negative relaxation at the 



problem depths - that is, we scale all corrections at those depths by 
a scale factor (typically 0.3). This procedure is often successful at 
accelerating convergence. 

(2) Ionizations fronts: These tend to show rapid changes over 
a few successive points. The position of the front oscillates, and 
there is often a strong coupling between the front location and the 
temperature. In these cases the most reliable technique to facilitate 
convergence is grid refinement in the neighborhood of the front. 
This is done semi-automatically in CMFGEN — a parameter file 
(which can be changed during the run) controls if, when and how 
the refinement is done. 



10 TESTING 

With a code as complicated as CMFGEN there is no easy way 
to fully test the code. However numerous checks have been per- 
formed to test individual components. For example, solutions ob- 
tained with different transfer modules have been compared, and 
found to agree (in limiting cases such as small v/c). At depth (e.g., 
r > 10), the solution to the gray and non-gray transfer equations 
also agree. Light curves of Type II SN computed with CMFGEN 
show reasonable a greement to those computed using th e hydrody- 
namics code V 1 D ( ILivnjl 19931 : [Pessart & Hillieill2010l) . A strong 
check, although not definitive, is also provided by our modeling 
of SN 1987 A. The time-sequence of models for SN 1987 A, be- 
gun at 0.3 d, shows a remarkabl e agreement to observations of SN 
1987A l lDessart&Hilliej2010l) — remarkable since there were no 
adjustable parameters in the modeling. 

Additional checks are provided by auxiliary files which list 
and check various processes and rates such as the ionization bal- 
ance of all species and the electron heating/cooling balance modi- 
fied for adv ection, adiabat ic cooling and nuclear energy deposition. 
As noted bv lHillieil ( l2003l) . the later equation, in the absence of SLs, 
is a linear combination of the radiative equilibrium equation and 
rate equations. In the presence of SLs this equality no longer holds, 
and thus the electron energy balance equation provides a check on 
possible problems with SL assignments. 



11 COMPUTATION OF OBSERVED SPECTRA 

The computation of the observed spectrum can be done in three 
ways: 

(i) Compute the observed spectrum using a Lorentz transforma- 
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tion of the comoving-frame boundary intensities computed by CM- 
FGEN (Section[TTTT)- 

(ii) Compute the observed spectrum using the comoving-frame 
boundary intensities computed using a separate calculation with 
CMF_FLUX. The advantage of this technique over the CMFGEN 
calculation is that we can use a finer spatial and frequency grid 
to improve the accuracy of the computation, we include all spec- 
tra l lines for the adopted mod el atoms, and we use the technique 
of lRvbicki&Hummej h994h to compute the electron-scattering 
source function. In this technique the redistribution of photons in 
frequency space due to the thermal motions of the electrons is ex- 
plicitly treated. The redistribution by the bulk motion of the gas is 
automatically allowed for by performing the radiative transfer cal- 
culations in the comoving frame. In the computation of the level 
populations (i.e., in CMFGEN) we generally treat the scattering co- 
herently in the comoving-frame - that is, we only allow for redistri- 
bution caused by the bulk motion of the gas. With the coherent scat- 
tering approximation the radiation is not explicitly coupled to that 
at other frequencies by a complex weighting function, and hence 
convergence is more stable. 

(iii) Using the emissivities and opacities computed in the co- 
moving frame, compute the observed spectrum using an observer's 
frame calculation. This calculation has a higher accuracy than the 
calculation done in the comoving frame. This occurs since, along a 
given ray, we are solving a differential equation in a single variable 
(z) whereas in the comoving frame we are solving a partial differ- 
ential equation in two variables (z,!/) (Section [l 1.2t . The differ- 
ence in accuracy is particularly noticeable for O stars - their narrow 
photospheric features become broadened as they are propagated (in 
both space and frequency) to the outer boundary. 

Ideally these techniques should give identical answers but in 
practice differences arise because the different techniques have 
different accuracies (as they use different numerical techniques) 
(Figs.m&H. 

At present we also ignore time dependence when computing 
observed spectra. This will be rectified in the future, but, for the 
same reasons outlined above (^, it is unlikely to have a large 
effect on observed spectra, except at the earliest times for H-rich 
core-collapse SNe when J and H are vaiying rapidly. For Type I 
SNe, the effects will be smaller since by one day the SN have al- 
ready cooled, and their spectra are slowly evolving. Note that time- 
dependence is taken into account in the calculation of the moments, 
and observer's frame spectra computed with these moments do dif- 
fer from those computed where the time-dependence of the mo- 
ments is ignored (Fig.|5}- The computation of the observer's frame 
spectrum is outlined in Section [Tl.2| 



11.1 Comoving Frame computation of tlie observed 
spectrum 

As we integrate from blue to red we store Io{t, Rmax, fJ^oji^o) 
at the outer boundary of the model. After the completion of 
each comoving frequency, we check whether we can compute 
/s(t, -Rmax, /is, i^s) for all Us for the next observer's frame fre- 
quency. If so, we compute Js(-Rmax, Mi fro™ hit, i?max, Moi 2^o) 

using the transformation from comoving to observer's frame 
(Ea.l42t and linear interpolation in frequency. Because of the way 
we constructed the ray grid there is a one-to-one correspondence 
between /io and /is, and hence no inteipolation in /i is necessary. 
The observed flux, Fs (t, i^s), at distance d is then computed using 



Fs{t,Ps) = -ip J pIs{t-d/c,Ryn!i^,p,Vs)dp. (37) 

This integration is performed using numerical quadrature with /is 
as the integration variable (note that pdp = -Rniax /J-sdus)- 



11.2 Observer's frame calculation of the observed spectrum 

To compute observed spectra we use a modified form of CMF_FLUX 
which has been described by^usche & Hillier (2005j. The calcu- 
lation of the observed spectrum proceeds as follows: 

(i) We solve the transfer equation in the comoving-frame us- 
ing the assumption of coherent electron scattering in the comoving 
frame. This allows us to compute the mean intensity in the comov- 
ing frame. Along each ray we insert extra grid points - typically we 
have a grid point every 0.5 Doppler velocities. The finer frequency 
grid adopted in CMF_FLUX generally has only a minor influence on 
the computed spectra for Typ e II SNe. 

(ii) Using the technique of iRvbicki & Hummed ( 1 19941) we com- 
pute the full, non-coherent and frequency dependent, electron- 
scattering source function. 

(iii) We resolve the transfer equation in the comoving frame, but 
this time we assume incoherent electron scattering, and hence use 
the electron-scattering source function previously computed. 

(iv) We repeat steps (ii) & (iii) until convergence is obtained. 
This generally only requires a few iterations although in model- 
ing the Type Iln SN 1994W a much larger number of iterations was 
needed. This was necessary because of the much higher line optical 
dept hs, at a given electro n density, induced by the lower SN veloc- 
ities toessart et ai]|2009l) . In most SNe, frequency redistribution to 
the red, associated with the bulk motion of the gas, dominates and 
the non-coherence of electron scattering due to thermal motions 
is subdominant. In SNe where the spectrum formation region ex- 
pands slowly, this non-coherent scattering is key and requires care- 
ful computation to yield accurate line-profile shapes. 

(v) We map the comoving frame emissivities and opacities 
into the obser ver's frame using the following relations (see, e.g., 
lMihalaslll978h : 



/3 = v/c (38) 

7 = (39) 

Xs(r,/is,!^s) = {uo/us)xoir,iyo) (40) 

r]sir,fis,Us) = [vs/voY ■no{r,Uo) (41) 

Isir,fis,Us) = {i^s/foflo{r,^o,yo) (42) 

Vs = Vo-yil + ^-oP) (43) 

Vo = Vsl{l - fJ.sP) (44) 



In the above we use "o" to denote the comoving frame and "s" to 
denote the observer's (static/rest) frame. This mapping provides the 
emissivities and opacities necessary to solve the transfer equation 
in the observer's frame using the usual (p, z) coordinate system. 
Along each ray (i.e., for each impact parameter p) we insert extra 
grid points so that the rapid variation of opacity and emissivity is 
adequately sampled - typically we have a grid point every 0.25 
Doppler velocities. 

(vi) We compute the observed flux using Eg. 1371 
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Figure 3. Comparison of the predicted observed spectrum computed using an observer's frame calculation (red), CMF_FLUX (blue), and CMFGEN (green) for 
a SN1987A like model at an age of 7.9 days. The main characteristics of the spectra are in very good agreement but the observer's frame calculation is of 
higher accuracy, and suffers less smoothing than spectra computed by mapping comoving-frame intensities, computed with CMF_FLUX and CMFGEN, into the 
observer's frame. 
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Figure 4. As for Fig.|3]but showing the logarithm of the flux over a broader wavelength range 
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12 GRAY TRANSFER 



we obtain 



To improve the speed of convergence it is desirable to have a good 
estimate of the temperature at each time step. In the inner, optically- 
thick regions, this is provided by solving the time-dependent gray 1 D{r J) ^ 1 d{r H) _ ^ — ^ J dv 
transfer problem. The solution of the time-dependent gray transfer Dt dr Jq 
also provides a check on the solution of the full-moment equations. f°° 



— XP {Si' - Jv)dv 
Jo 

Integrating Eqs.|5]and|6]over frequency, and regrouping terms, ~ xp{S — J) (45) 
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Figure 5. Comparison of the predicted observed spectrum computed allowing for time dependence in the moment equations (red) with the solution for the 
moments computed using the fully relativistic solution but ignoring time-dependence (blue). Both spectra were computed in the observer's frame. 



and 



1 D{r-^H) 1 d{r'^K) K - J 



Dt 



+ 
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+ 



-XR I 

'0 



-XrH 



(46) 



where XR is the Rosseland-mean opacity, and xp is the Planck- 
mean opacity. Assuming LTE holds at depth we have = 
Using equation|27]with|45]we obtain 



1 D(r*J) 1 d{r^H) 



Dt 



+ 



dr 



. I f^dccay 



De P Dp\ 
(47) 

which follows from Eq.|27] After solving for J we need to relate J 
to B. Using Eq.|27] and since S = B,we have 



X^B^dv = 



XvJi^dv 



+ 



1 

47r 



De P Dp 



Dt 



and thus 



B = aT* 



Atixr 



fide 



De P Dp 



(48) 



(49) 



where we have assumed Ji, ~ . 

The solution of Eq. |47]is complicated by the presence of the 
pDe/Dt — PDln p/ Dt (= W) term which depends on the un- 
known temperature structure. We handle the term by iteration — 
that is, we compute the temperature structure for a given W by 
solving Eqs. 147 1 and 1491 update W, and then resolve Eqs. 147 1 and 
|49l Because of stability issues in the outer regions, we found it 
necessary to limit the changes in W to achieve convergence. To de- 
termine the change in e we assume that the departure coefficients 
are constant - in practice this is equivalent to assuming UTE in the 
inner regions where the gray approximation is valid. 



Equations|46l|47] and|49|represent 3 equations in 4 unknowns 
(J, H, K, T). As for the frequency-dependent transfer equations we 
eliminate K using the Eddington relation K = fj where / is ob- 
tained by solving the time-independent transfer equation, and for 
simplicity we neglect the curvature of the characteristics. 

Since the Rosseland-mean opacity also depends on tempera- 
ture, we perform several iterations of the gray-solution with up- 
dated Rosseland-mean opacities. As the gray solution only provides 
a starting solution, high accuracy is unnecessary, and practical ex- 
perience shows that the above procedure provides a very acceptable 
starting solution. As the gray solution does not hold everywhere, we 
only apply the gray-temperature structure to regions with an opti- 
cal depth above Tmin- Below Tmin we use T from the previous time 
step. Around Tmin we use a simple interpolation procedure so as 
to switch smoothly between the two approximations. The choice 
of Tmin is somcwhat problematical. In O and WR stars we found 
that it was reasonable to adopt Tmin ~ 1 however for Type II SNe 
we found that Tmin ~ 5 (or 10) provided better starting solutions. 
A comparison between the actual model temperature structure and 
the gray temperature structure is show in Fig. |6] 



12.1 Global Energy Constraint 

A global energy constraint can be obtained by integrating Eq. l l47t 
over volume. Multiplying Eq. |47]by and integrating over r we 
obtain 



+ 



rH{r) 

J_ D{r^J) 
cr^ 



Dt 



dr 



De P Dp 



Dt 



(50) 



In a static atmosphere, with v = this simply reduces to 
r^iif =constant; and this is often used as the energy constraint in 
stellar atmosphere calculations at large r (r > 1). This constraint 
has real physical meaning — the flux at the outer boundary is set 
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Figure 6. Comparison of the profiles for the initial gray temperature Tq 
and the final model temperature Tp. At depth the agreement is excel- 
lent. The maximum difference between the two temperatures is about 
10%. At V =10000kms"i the Rosseland optical depth is ~ 12, while 
at 30000kms-l it is ~0.00035. Corrections for v >30000kms-l are 
small because this model has an enforced minimum temperature. 



by the flux imposed at the inner boundary. However, for time- 
dependent SN calculations at early times, the constraint equation 
is, in some sense, merely a mathematical constraint. Because of the 
finite speed of light coupled with the effects of diffusion, the outer 
flux is independent of the inner atmosphere at the current time step 
when the envelope has a "large" optical depth. Nevertheless, Eg. 1501 
does allow a consistency check, and in complex codes consistency 
checks are always useful. In the outer envelope, and as the envelope 
thins, and in Type la SN with their much lower envelope masses, 
the constraint represented by Eq. |50] becomes much more phys- 
ically meaningful. Of course, depending on the balance between 
adiabatic cooling and other heating and cooling processes, errors 
in the temperature of the inner envelope may manifest themselves 
at a later time. 

In our original formulation, we integrated from rmin to r. 
However, we found such a formulation was not so useful since nu- 
merical errors (at large optical depths) tend to limit its usefulness. 
These numerical eiTors arise from the discretization {J^ is defined 
on the grid while Hi, is defined at the midpoints), from the eval- 
uation of the integrals, from the need for very high convergence, 
and from the technique used to difference the moment equations. 
For example, when we difference the zeroth moment equation we 
effectively replace 



at each depth d by 



/ Xi'{Sv - J„)du 
Jo 



< Xv > {Sv — Jv)dv 



where < Xi' > is now averaged over several depthfl By integrat- 
ing inwards, we can more directly compare the constancy of the 
right-hand-side of equation |50] with the quantity of interest — the 
flux (or luminosity) at the outer boundary. 



^ The averaging in x occurs because of the definition of the optical depth 
increments. 



In Fig. m and Fig. [H we illustrate the global energy constraint 
as a function of depth for a Type la SN model, and model for SN 
1987A. In the outer regions, constancy of the RHS is excellent, 
although small departures in the inner regions can be seen. These 
departures often arise from sharp changes in the material properties 
- e.g., the H ionization front in Type II SN, or from sharp changes 
in composition. By improving the resolution across these features 
we improve the accuracy of the global energy constraint. Through 
extensive testing we found that, in general, improving the accuracy 
with which the global energy constraint was satisfied did not sig- 
nificantly change predicted spectra, especially at earlier epochs. 



13 FUTURE WORK 

We have described the procedure we use to undertake time- 
dependent modeling of homologous SN ejecta. With our work, 
we solve self-consistently the time-dependent radiative-transfer 
equation and its 0th and 1st moments, the time-dependent energy 
balance equation, and the time-dependent kinetic rate equations. 
Our technique has been successfully used to model the Type II- 
pecuhar SN I987A jPessart & Hillieill2010l) . Type I l-Plateau SNe 
I Dessart &Hillieill201ll) , and Type Ilb/Ib/Ic SNe jPessart et af 



20111). We are also using it to model Type la SNe i Dessart et al, 
20I2al) and extend our modelling of Type Ib/Ic SNe jDessart et al 
20I2bl) . 

With recent advances in hydrodynamic simulations there is in- 
creasing evidence that SNe are not spherically symmetric or homo- 
geneous. There is also observational evidence for non-symmetric 
ejecta and for mixing in SNe. Ultimately 3D radiative transfer 
models may be needed to accurately model SNe ejecta. How- 
ever, for a variety or reasons, ID models still have a very impor- 
tant role. As our own work has shown, ID models can be used 
to identify important fundamental physics such as the importance 
of ti me dependence in the rate equations for modeling Type II 
SN jDessart & Hillieill2008l) or the prediction that He I emission 
lines can be excited in ea rly spectra of lb S N without the need 
for non-thermal processes dDessart et al.ll20I H) . They also can pro- 
vide fundamental insights into the abundance and physical struc- 
ture of ejecta - iirformation which is missing for most SN ejecta. 
Discrepancies between ID models and observations can also pro- 
vide crucial insights into the need for missing physics, and for fea- 
tures which can only be explained by 3D effects. Since hydrody- 
namical modeling cannot yet accurately predict the 3D structure of 
ejecta, 3D modeling is necessarily rich in free-parameters which 
can limit crucial insights. Importantly, ID models also provide a 
testbed against which more complex models can be developed. 

The present models are restricted to the case of SN ejecta in 
homologous expansion. This is of little consequence for Type la 
SNe, but restricts modeling for Type Il-peculiar objects like SN 
1987 A to t > 1 d, and to even later times {t > lOd) for "generic" 
Type II SNe of the plateau type. In order to overcome this limita- 
tion, work on a fuU-relativistic time dependent solver is in progress. 

Another issue which we are addressing is the influence of non- 
thermal ionization and excitation. Such processes are believed to 
be i mportant fo r explaining the appearance of He I lines in Type lb 
SNe l iLucvlfl991, ; . Dessart et al. 201 1, 2012bi) , and for the produc- 
tion of H I lines in the nebular phase of Type II S Ne, and in 1987A 
jXu & McCra"vlll99lLlKozma & Franssonlll998j|bh . 

In our work we treat the whole ejecta. For Type I SNe, with 
their relatively low ejecta masses, this is a reasonable approach. 
However, for Type II SNe, with their large ejecta masses, our ap- 
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Figure 7. Illlustradon of the global energy constraint for a Type la SN model with local energy deposition at an age of 60 days. In the left plot we see the 

depth-dependent luminosity (red; solid), the "conserved luminosity" (blue; ), and the contiibutions to this conserved quantity by energy from radioactive 

decay (green; . — . — .— ), from the Dr"^,]/ Dt term (purple; ...). The gas term is neghgible and is not shown. In the right plot we see the percentage error 
in the conserved luminosity (red; solid), and the total correction to the luminosity as a function of depth, normalized by the outer boundary comoving-frame 
luminosity (blue; ). The model was computed using 114 depth points. 
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Figure 8. Illlustration of the global energy constraint for a SN 1987A-like model with local energy deposition at an age of 37 days. In the left plot we see the 

depth-dependent luminosity (red; solid), the "conserved luminosity" (blue; ), and the contributions to this conserved quantity by energy from radioactive 

decay (green; . — . — ), Dr^J / Dt term (purple; ...), and the gas term (pink;—...). In the right plot we see the percentage error in the conserved luminosity 

(red; solid), and the total correction to the luminosity as a function of depth, normalized by the outer boundary comoving-frame luminosity (blue; ). 

The model was computed using 126 depth points. 



proach may be an overkill. The inner layers of these models can 
be adequately described using a time-dependent LTE gray model. 
A superior approach might be to combine a gray model with CM- 
FGEN modeling for optical depths r < rmax (where Tmax may be 
of the order of 100), and use the gray solution to provide a lower 
boundary condition for CMFGEN. This would reduce the number of 
depth points in the CMFGEN models, and potentially the computa- 
tional effort (since the model computation time scales as the square 
of the number of grid points). 

With these tools we will be able to address many important 
questions related to SNe, their progenitors, and their mechanisms 
of explosion. For example, do stellar-evolution and hydrodynami- 
cal models of SNe agree with observations? What progenitors give 
rise to what SN types? Do SN Ilb/Ib/Ic arise from binary-star evolu- 
tion? Does the composition of the outer unbumt ejecta correspond 
to predictions of stellar evolutionary models and stellar- atmosphere 
calculations of their progenitor class? Do we really understand the 
homogeneity of Type la SNe and the origin of the width-luminosity 
relation? Is there some signature in the spectra of Type la SNe that 



may distinguish their origin from a single- or double-degenerate 
scenario? 



APPENDIX: 7-ray MONTE-CARLO TRANSPORT FOR SN 
EJECTA 

We have two possibilities for the treatment of the decay energy 
in CMFGEN. We may deposit the entire decay energy locally and 
neglect any 7-ray transport, as described in Section |5] This ap- 
proximation is very good for hundreds of days in Type II SNe be- 
cause unstable isotopes are produced explosively within the slower- 
expanding higher-density innermost layers of the massive ejecta. 
This is partially supported by the nebular-phase observations of SN 
1987A , which show the luminosity fading expected for full 7 -ray 
trapping for up to ^^^SOOd after explosion jAmett et al.lll989l) . On 
the other hand 7-rays were detected as early as day 200, indicat- 
ing that some 7-rays were escaping at that epoch, and by inference, 
that some non-local deposition was also occurring at that epoch. 
This, and the subsequent evolution of the 7-rays and X-rays from 
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SN 1987A, provide strong constraints on the r nixing that has oc- 
curre d (e.g., Kumagai et al. ' 1989*: ' Xrnett & FulfigM IFu & AmettI 
[1989. : Burrows & van Riper 1995,). 

However, in the lower mass ejecta that characterize Type I 
SNe, either of thermonuclear or core-collapse origin, non-local 
energy deposition and even 7-ray esc ape can occur on a much 
shorter time scale dOessart et al.ll201lh . In Type la SNe, the de- 
flagration/detonation produces ^®Ni at much larger velocities than 
in core-collapse SNe, in layers that can be well above the eject a 
base depending on the details of the explosion ([Kho^io3[l99l|). 
This p roperty leads to 7-ray escape as early as the peak of the light 
curve jHoeflich et al.|[l99a) . Hence, to model such SN ejecta, one 
really needs to solve for the transport of 7-rays produced in the 
decay of unstable isotopes. 

Since our Monte Carlo code is to provide an energy-deposition 
distribution for CMFGEN runs, it uses the same model atom (i.e., 
species and ions), the same spatial grid, and assumes homolo- 
gous expansion. The code treats only two-step decay chains but 
as many as desired, as in CMFGEN. The main focus for now 
though is on the decay chain associated with ^®Ni. The decay 
lifetimes, 7-ray line energies (and probabilities), and positron en- 
ergies are all taken from http ://w w w . nndc .bnl . gov/chart/ [ them- 
selves based on the work o f IHuo et all l fT987h . For ^"Ni and '^'^Co 
decays, these are in close ag r eement with the v alues given by 
lAmbwani & Sutherland! ( Il988l) ; iNadvozhinI ( Il994l) . The total de- 
cay energy for each channel, obtained by summing over all 7-ray 
lines and including positron energy, is used in CMFGEN when as- 
suming local energy deposition. For non-local energy deposition, 
the Monte Carlo code is used instead to treat all individual contri- 
butions. This ensures energy consistency whether we assume local 
or non-local deposition. 

At present, we do not use the Monte Carlo procedure to select 
7-ray lines but instead treat all important 7-ray photons associated 
with a given decay. This is more time consuming but produces bet- 
ter statistics for weaker and/or improbable 7-ray lines. In practice, 
to reduce the burden on the number of 7-ray lines we follow, we 
treat only 7-ray lines with a probability higher than 1%. This cor- 
responds to 6 7-ray lines for the decay of an ^^Ni nucleus and to 
15 for ^®Co. For ^''Co decay, the missing 7-ray lines cause a <3% 
deficit in the decay energy, which we compensate by scaling all 
^^Co 7-ray line energies accordingly. 

The Mo nte Carlo procedure we follow is analogous to that 
described by lAmbwani & Sutherland' ("19881) while the numerical 
approach closely follows Hillier (1991). Although our SN mod- 
els assume spherical symmetry, our calculations are performed on 
a 3D spatial grid. In the Monte Carlo technique, one has first to 
build probability functions for all important processes. For a given 
event described by the function / of a variable x over the interval 
[Xmin, 2;inax], wc cau generate x from 

q = f fix)dx/ f fix)dx (51) 

^min ^min 

w here 17 is a ra ndom number uniformly distributed between and 
1 ( lHillieill99ll) . One can also bias the sampling of f{x) by a user- 
specified function b{x) so that x is drawn from 

q= b{x)f{x)dx/ / h{x)f{x)dx (52) 

^min ^miii 

and that decay is then given a weight l/b[x). 

We thus build probability distributions to describe the total 
number of decays, the relative number associated with each 2-step 
decay chain considered, the relative number associated with each 



nucleus in each chain, and finally the radial distributions of such 
decays. 

For the simulation of 7-ray spectra, it is advantageous to bias 
the sampling of decaying nuclei to favor emission from the outer 
ejecta regions and thus enhance the statistics (reduce the noise) of 
the emergent 7-ray spectrum. This can help at early times when the 
bulk of unstable isotopes are located at large optical depths so that a 
very small fraction of 7-rays manage to escape. When biasing 7-ray 
emission for escape, we use the bias function h(r) = exp(— r(r-)) 
where r(r) = J^°° K{r')p{r')dr' . For simpli city, we use a fixed 
mass- absorption coefficient k = 0.03 cm^ g~^ l lKozma & Franssonl 
I1992I : this value does not need to be accurate). The associated en- 
ergies for this decay at r are then weighted by l/fe(r). In practice, 
this option is rarely used. Indeed, we are interested by the energy 
fraction that is deposited, not that which escapes. Furthermore, the 
eventual detection of 7-rays from SNe will most likely occur for an 
extragalactic event when the bulk of the released energy escapes as 
the ejecta becomes transparent to 7-rays — biasing emission is no 
longer needed at such times. 

Using a series of random numbers g, we select the decay 
chain, the nucleus that decays, and the decay location in spheri- 
cal coordinates (r, 6, cj)). r is selected according to Fan. 15 II with 
X = r and f{x) — Anr'^r} where 77 is the emissivity. The angles are 
selected using 

cos 9 = 2q ~ 1 , and (j) = ""(2? — 1) , 

and a similar procedure is used to determine the original orientation 
of 7-ray emission. 

We then proceed to transport all 7-ray photons associated 
with that decay (note that we neglect the expansion of the ejecta 
over the life of such 7-rays — the expanding ejecta appears as 
frozen to those 7-rays). If the decay also produces positrons, these 
are deposited locally as heat. The location of the next scatter- 
ing/absorption is conditioned by the opacity. We consider the two 
processes of Compton scattering and photoelectric absorption (we 
neglect the opacity as sociated with pair p roduction) and adopt the 
analytical formulae of lKasen et alj ^20061) . Given an orientation we 
compute the associated total optical depth along that ray and direc- 
tion out to the outer ejecta boundary. We assume a hollow core for 
rays that cross the inner-boundary of the ejecta. Using another ran- 
dom number q, the location of the next scattering/absorption is at an 
optical depth r given by g = 1 — exp(— r) or r = — log(g — 1). In 
calculating these ray optical depths, we account for Doppler effects 
(redshifts) associated with expansion and include both scattering 
and absorption opacities. If r is greater than the maximum optical 
depth along the ray in the selected direction, the 7-ray photon es- 
capes. Otherwise, we step along the ray to find the new location, 
change the photon energy to the scatterer's frame, and generate an- 
other random number to determine whether the photon is absorbed 
or scattered at that location. If the latter case is drawn, a new direc- 
tion of scattering is sampled from the differential cross section for 
Compton scattering and the new energy of the photon is computed 
(the energy loss in the scattering is then deposited locally). 

In our calculations, we typically use 100,000 decays, which 
coTOsponds to about ten times as many 7-ray photons depending 
on the post-explosion time (i.e., whether it is primarily ^''Ni or 
^''Co nuclei that decay). The agreement between the local (com- 
puted analytically) and non-local (computed numerically with the 
Monte Carlo code) energy deposition profiles at very early times 
agree to within a fraction of a percent in most regions. The two 
can depart in regions where the original ^''Ni mass fraction is very 
low since too few decays will be triggered in those regions. Since 
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Figure 9. Bolometric 7-ray light-curve fro m >10d until 1500 d afte r ex- 
plosion for a SN la model (the model from lKasen & WooslevI )2007l) with 
0.49 Mn of ^"^N i initia lly) and a SN II-P model (the model sl5el2 from 
iDessart & HillieJ i201ll) . endowed with 0.087 Mq of s^Ni initially). We 
show the evolution of the total decay energy (solid), the total decay en- 
ergy deposited in the ejecta (dash-dotted), and the escaping 7-ray energy 
(dashed). The orange line traces the radioactive-decay power of 1 Mq of 
^^Ni produced initially in an explosion. 



we are mostly concerned with the bulk of the energy that is effec- 
tively deposited (and in appreciable amounts to influence the gas 
properties), we are confident that the (originally) ^^Ni-rich regions 
are well sampled and the actual distribution profile is accurately 
determined. 

We illustrate some interesting observable properties pro- 
duced by such transport calculations. In Fig. |9] we illustrate 
the evolution of the total, deposited, and escaping decay energy 
fro m "rep resentative" SN l a (the model with 0.49 Mq of ^®Ni 
of iKasen & Wogsle^] l2007h and II-P (the model sI5el2 from 
lOessart & HillieJ 201 1 ) models. In the Type II-P model, the larger 
mass, the slower expansion and the deeper location of ^^Ni pro- 
duced by the explosion lead to a much delayed emergence of 7-rays 
at a few hundred days after explosion. In contrast, the SN la model 
lets 7-rays start to escape as early as early as two weeks after ex- 
plosion. In each case, we illustrate the 7-ray spectrum one would 
observe at an epoch when 7-ray escape becomes strong. Because 
of the relatively short half-life of ^''Ni nuclei, the 7-ray lines stem 



primarily from ' Co nuclei decay. 

We do not treat time-dependence in the Monte-Carlo code. 
This is a reasonable approximation since the light-travel time is 
small compared with the flo w time, and because gamma-rays do 
not diffuse. As discussed bv lSwartz et al.| Il995b a 7-ray deposits 
half if its energy, on average, per scattering, and hence, after a few 
scatterings most of the gamma-ray energy has been deposited. 
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